/*
 * Assemble.cpp
 *
 *  Created on: 1 Aug 2011
 *      Author: Allan
 */

#include "Assemble.h"

// C++ includes
#include <algorithm>
#include <set>

// Boost includes
#include <boost/foreach.hpp>

//GeoReact includes
#include "Database.h"
#include "EquilibriumConstantData.h"
#include "SpeciesMap.h"

#include "AqueousSolution.h"
#include "GaseousSolution.h"
#include "MineralSolution.h"

#include "Phase.h"
#include "Multiphase.h"

#include "Reaction.h"
#include "ReactionKinetics.h"
#include "ReactionSystem.h"

#include "Kinetics/MineralKinetics/MineralKineticsDatabase.h"

const Multiphase AssembleMultiphase(const vector<string>& allSpecies)
{
	vector<Phase> phases;
	
	vector<string> aqueousSpecies = FilterSpecies(allSpecies, {"aqueous"});
	vector<string> gaseousSpecies = FilterSpecies(allSpecies, {"gaseous"});
	vector<string> mineralSpecies = FilterSpecies(allSpecies, {"mineral"});
	
	if(!aqueousSpecies.empty())
		phases.push_back( Phase(AqueousSolution(aqueousSpecies)) );
	
	if(!gaseousSpecies.empty())
		phases.push_back( Phase(GaseousSolution(gaseousSpecies)) );
	
	BOOST_FOREACH(const string& mineral, mineralSpecies)
		phases.push_back( Phase(MineralSolution({mineral})) );
	
	return Multiphase(phases);
}

///====================================================================================================///
/// Methods for assemblage of reactions
///====================================================================================================///
const Reaction AuxiliarAssembleReaction(const ReactionInfo& reactionInfo)
{
	// Create the reaction object
	Reaction reaction;
	
	// Set the attributes of the reaction and then register its index related data
	reaction.SetEquation(reactionInfo.reactant, reactionInfo.products, reactionInfo.stoichiometries);
	reaction.SetEquilibriumConstantData( EquilibriumConstantData(reactionInfo.logKi, reactionInfo.Ti) );
	
	// Return the reaction object
	return reaction;
}

const ReactionSystem AssembleReactionSystem(const vector<string>& species)
{
	vector<Reaction> reactions; reactions.reserve(species.size());
	
	BOOST_FOREACH(const string& s, species)
		reactions.push_back( AssembleReaction(s) );
	
	return ReactionSystem(reactions);
}

const Reaction AssembleReaction(const string& species)
{
	if(Database::GetInstance().GetSpeciesInfo(species).HasTag("basis"))
		return AssembleBasisReaction(species);
	
	else if(Database::GetInstance().GetSpeciesInfo(species).HasTag("mineral"))
		return AssembleMineralReaction(species);
	
	else 
		return AssemblePlainReaction(species);
}

const Reaction AssemblePlainReaction(const string& species)
{
	// Get the information data of the reaction whose reactant is the given species
	const ReactionInfo reactionInfo = Database::GetInstance().GetReactionInfo(species);
	
	return AuxiliarAssembleReaction(reactionInfo);
}

const Reaction AssembleBasisReaction(const string& species)
{
	ReactionInfo reactionInfo;
	
	reactionInfo.reactant        = species;
	reactionInfo.products        = {species};
	reactionInfo.stoichiometries = {1.0};
	reactionInfo.logKi           = {0.0};  // K = 1
	reactionInfo.Ti              = {25.0}; // T = 25 degC

	return AuxiliarAssembleReaction(reactionInfo);
}

const Reaction AssembleMineralReaction(const string& mineral)
{
	// Create the reaction object for the given mineral
	Reaction mineralReaction = AssemblePlainReaction(mineral);
	
	// Look up at the database of mineral kinetics for the given mineral
	optional<MineralKineticsInfo> mineralKineticsInfo = MineralKineticsDatabase::GetInstance().FindMineralKineticsInfo(mineral);
	
	// If the kinetics information of the minaral was found, create a mineral kinetics intance and set the mineral reaction with it
	if(mineralKineticsInfo)
	{
		// Create the mineral kinetics object of the given mineral
		MineralKinetics mineralKinetics;

		// Set the mineral kinetics object with the found information
		mineralKinetics.SetReaction(mineralReaction);
		mineralKinetics.SetSpecificSurfaceArea(mineralKineticsInfo->specificSurfaceArea);
		mineralKinetics.SetMechanisms(mineralKineticsInfo->mechanisms);
		
		// Set the kinetics of the mineral reaction
		mineralReaction.SetReactionKinetics( Reaction::ReactionKineticsPtr(new MineralKinetics(mineralKinetics)) );
	}
	
	return mineralReaction;
}

///====================================================================================================///
/// Methods for assemblage of species
///====================================================================================================///
const vector<string> AssembleSpecies(const vector<string>& knownSpecies)
{
	return AssembleSpecies(knownSpecies, {});
}

const vector<string> AssembleSpecies(const vector<string>& knownSpecies, const vector<string>& tagsSecondarySpecies)
{
	const vector<string> basisSpecies = AssembleBasisSpecies(knownSpecies);
	const vector<string> secondarySpecies = AssembleSecondarySpecies(knownSpecies);
	
	vector<string> filteredSecondarySpecies = secondarySpecies;
	
	// Filter the secondary species that contains the given tags (if there is any)
	if(!tagsSecondarySpecies.empty())
		filteredSecondarySpecies = FilterSpecies(secondarySpecies, tagsSecondarySpecies);
	
	set<string> allSpecies;
	
	allSpecies.insert(basisSpecies.begin(), basisSpecies.end());         
	allSpecies.insert(filteredSecondarySpecies.begin(), filteredSecondarySpecies.end());
	allSpecies.insert(knownSpecies.begin(), knownSpecies.end()); // Make sure the secondary species in (knownSpecies) are inserted as well
	
	return vector<string>(allSpecies.begin(), allSpecies.end());
}

const vector<string> AssembleBasisSpecies(const vector<string>& knownSpecies)
{
	// Get an instance of the database
	const Database& db = Database::GetInstance();
	
	// The set containing only unique basis species
	set<string> basisSpecies;
	
	// Make sure the hydrogen ion H[+] is present in the set of basis species in case H2O is known to exist
	if(std::find(knownSpecies.begin(), knownSpecies.end(), "H2O") != knownSpecies.end())
		basisSpecies.insert("H[+]");
	
	// Loop over all the known species
	BOOST_FOREACH(const string& species, knownSpecies)
	{
		// Check if the current known species is a basis species in the database
		if(db.GetSpeciesInfo(species).HasTag("basis"))
			basisSpecies.insert(species);
		else
		{
			// Get the info data of the reaction whose reactant is the current known species
			const ReactionInfo reactionInfo = db.GetReactionInfo(species);
		
			// Add all the products in the reaction to the set of basis species being created (all these products are basis species in the database)
			basisSpecies.insert(reactionInfo.products.begin(), reactionInfo.products.end());
		}
	}
	
	// Convert the set of basis species to a vector container and return it
	return vector<string>(basisSpecies.begin(), basisSpecies.end());
}

const vector<string> AssembleSecondarySpecies(const vector<string>& knownSpecies)
{
	const vector<string> basisSpecies = AssembleBasisSpecies(knownSpecies);
	
	// Get an instance of the database
	const Database& db = Database::GetInstance();
	
	// Get a reference to the list of reactions in the database
	const Database::ReactionInfoList& reactionInfoList = db.GetReactionInfoList();
	
	// The vector of secondary species that can be built from the set of basis species
	vector<string> secondarySpecies;
	
	// Loop over all the reactions in the database
	BOOST_FOREACH(const ReactionInfo& r, reactionInfoList)
	{
		// An auxiliary boolean variable to determine if all the products of the current reaction belong to the set of basis species
		bool belong = true;
		
		//  Check if the products of the current reaction belong to the set of basis species 
		BOOST_FOREACH(const string& product, r.products)
			if(std::find(basisSpecies.begin(), basisSpecies.end(), product) == basisSpecies.end())
				belong = false;
		
		// Insert the reactant (a secondary species) of the current reaction into the vector of secondary species
		if(belong)
			secondarySpecies.push_back(r.reactant);
	}
	
	return secondarySpecies;
}

const vector<string> AssembleSecondaryAqueousSpecies(const vector<string>& knownSpecies)
{
	const vector<string> secondarySpecies = AssembleSecondarySpecies(knownSpecies);
	
	const vector<string> aqueousSecondarySpecies = FilterSpecies(secondarySpecies, {"aqueous"});
	
	return aqueousSecondarySpecies;
}

const vector<string> FilterSpecies(const vector<string>& species, const vector<string>& tags)
{
	// The subset of the given species that have any of the given tags
	vector<string> newSpecies;
	
	// Loop over all the given species 
	BOOST_FOREACH(const string& s, species)
		if(Database::GetInstance().GetSpeciesInfo(s).HasAnyTag(tags))
			newSpecies.push_back(s);
			
	return newSpecies;
}

const vector<string> RemoveSpecies(const vector<string>& species, const vector<string>& tags)
{
	// The subset of the given species that does not have any of the given tags
	vector<string> newSpecies;
	
	// Loop over all the given species 
	BOOST_FOREACH(const string& s, species)
		if(!Database::GetInstance().GetSpeciesInfo(s).HasAnyTag(tags))
			newSpecies.push_back(s);
			
	return newSpecies;
}
